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ABSTRACT 
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A finite difference method for the solution of symmetric posi- 
tive linear differential equations is developed. The method is 
applicable to any region with piecewise smooth boundaries. Methods 
for solution of the finite difference equations are discussed. The 
finite difference solutions are shown to converge at essentially 
the rate CKh 1 / 2 ) as h -*■ 0, h being the maximum distance between 
adjacent mesh points. 

An alternate finite difference method is given with the ad- 
vantage that the finite difference equations can be solved itera- 
tively. However, there are strong limitations on the mesh arrange- 
ments which can be used with this method. 

The Tricomi equation can be expressed in symmetric positive 
form. Admissible boundary conditions for any region with piece- 
wise smooth boundaries are given, with a wide choice of boundary 
conditions being possible. 

A Tricomi equation with a known analytical solution is solved 
numerically as an illustration of the numerical results which can 
be obtained. There is strong convergence to the analytical solu- 
tions, but pointwise divergence. Smoothing of the solution reduces 
this, though, and satisfactory numerical results are obtained. 
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INTRODUCTION 


In the theory of partial differential equations there is a 
fundamental distinction between those of elliptic, hyperbolic and 
parabolic type. Generally each type of equation has different 
requirements as to the boundary or initial data which must be 
specified to assure existence and uniqueness of solutions, and to 
be well posed. These requirements are usually well-known for an 
equation of any particular type. Further, many analytical and 
numerical techniques have been developed for solving the various 
types of partial differential equations, subject to the proper 
boundary conditions, including even many nonlinear cases. However, 
for equations of mixed type much less is known, and it is usually 
difficult to know even what the proper boundary conditions are. 

As a step toward overcoming this problem Friedrichs [l] has 
developed a theory of symmetric positive linear differential equa- 
tions independent of type. Chu [2] has shown that this theory can 
be used to derive finite difference solutions in two -dimens ions for 
rectangular regions, or more generally, by means of a transformation, 
for regions with four corners joined by smooth curves. In this 
paper a more general finite difference method for the solution of 
symmetric positive equations is presented. The only restriction on 
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the shape of the region is that the boundary be piecewise smooth. 
It is proven that the finite difference solution converges to the 
solution of the differential equation at essentially the rate 
0(bV 2 ) as h -* 0, h being the maximum distance between adjacent 
mesh points for a two-dimensional region. Also weak convergence 
to weak solutions is shown. 

An alternate finite difference method is given for the two- 
dimensional case with the advantage that the finite difference 
equation can be solved iteratively. However, there are strong 
limitations on the mesh arrangements which can be used with this 
method. 

As an example of the potential usefulness of the theory of 
symmetric positive equations, the Tricomi equation 


y^xx ■ v 

can be expressed in symmetric positive form. It is shown that 
suitable boundary conditions can always be determined, regardless 
of the shape of the region. The problem in a practical case is to 
determine an "admissible" boundary condition which corresponds to 
available boundary information. 

As an illustration of numerical results which can be obtained 
by the proposed finite difference scheme, a Tricomi equation with 
a known analytical solution is solved numerically. The results in- 
dicate that, although there is strong (i.e., L^) convergence of the 
finite difference solution to the analytical solution, there is 
pointwise divergence along the boundary. However, smoothing the 
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solution can eliminate this problem, and satisfactory numerical 
results are obtained, although rigorous mathematical Justification 
of the smoothing process is not given. 



CHAPTER I 


SYMMETRIC POSITIVE LINEAR DIFFERENTIAL EQUATIONS 
1.1 Basic Definitions 

Let ft be a bounded open set in the m-dimensional space of 
real numbers, R m . The boundary of ft will be denoted by dft, and 
its closure by ft. It is assumed that dft is piecewise smooth. 

A point in R m is denoted by x = (x^xg, . . . , x^j) and an 
r-dimensional vector valued function defined on ft is given by 
u = (u-^Ug, . . ., u,,). Also let a?-, a?, . . . , a m and G be 
given r X r matrix -valued functions and f = (f^fg, . . . , f r ) 
a given r dimensional vector-valued function, all defined on ft 
(at least). It is assumed that the aA are piecewise differen- 
tiable. For convenience, let a = (a?-, a?, . . ., a m ), so that we 
can use expressions such as 

m n 

V • ( a u) = J ( a *u) . (l.l) 

i=l 

With this notation we can write the identity 

(aM = ^^ u+ ^ al S 7 

i=l i=l i=l 

simply as 

V • (au) = (V • a) u + a ' Vu (l.2) 
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With this we can give the definitions for symmetric positive 
operators and admissible or semi -admissible boundary conditions 
which were introduced by Friedrichs [l] . 

Let K be the first order linear partial differential opera- 
tor defined by 

Ku = a • Vu + V • (au) + Gu (1.3) 

K is symmetric positive if each component, a^, of a is symmetric 
and the symmetric part, (G + G*)/2, of G is positive definite on 
on f i. 


For the purpose of giving suitable boundary conditions, a 
matrix, 3> is defined (a.e.) on bQ by 

3 = n • a (1.4) 

where n = (n^ng, . . ., r^) is defined to be the outer normal 
on bQ. 


The boundary condition Mu = 0 on dfi is semi -admissible 


if M = p - 3 , where p is any matrix with non-negative definite 
symmetric part, (p + p*)/2. If in addition, H(p - 3)®Tl(p + 3 ) = R r 
on the boundary, bQ, the boundary condition is termed admissible . 
(Tl(p - 3 ) is the null space of the matrix (p - 3 ).) 

The problem is to find a function u which satisfies 
Ku = f on G 

Mu = 0 on bQ 

where K is symmetric positive. 

It turns out that many of the usual partial differential equa- 


} 


(1.5) 


tions may be expressed in this symmetric positive form, with the 
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standard boundary conditions also expressed as an admissible bound- 
ary condition. This includes equations of both hyperbolic and ellip- 
tic type. However ; the greatest interest lies in the fact that the 
definitions are completely independent of type. An example of 
potentially great practical importance is the Tricomi equation 
which arises from the equations for transonic fluid flow. The 
Tricomi equation is of mixed type, i.e., it is hyperbolic in part 
of the region, elliptic in part, and is parabolic along the line 
between the two parts. 

The significance of the semi -admissible boundary condition 
is that this insures the uniqueness of a classical solution to 
a symmetric positive equation. On the other hand, the stronger, 
admissible boundary condition is required for existence. The 
existence of a classical solution is generally difficult to prove 
for any particular case, and depends on properties at corners of 
the region. However, it is very easy to prove existence (but not 
uniqueness'.) of weak solutions with only semi -admissible boundary 
conditions . 

1.2 Basic Identities and Inequalities 

Let 34 be the Hilbert space of all square integrable 
r-dimensional vector-valued functions defined on fl. The inner 
product is given by 


(u,v) = f 


U “ V 


( 1 . 6 ) 
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where 


u • 



and 

||u|| 2 = (u,u) 

A boundary inner product is defined by 


(1.7) 


(u,v 



U • V 


( 1 . 8 ) 


with the corresponding norm 

M|| = (u,u) B (1.9) 

We introduce now the adjoint operators K* and M*, which are 
defined by 

K*u - -a • Vu - V • (au) + G*u (l.lQ) 

M*u = (p* + 3)u (l.ll) 

The relation between K and M and their ad joints is given 
by Friedrichs "first identity." 

Lemma 1.1 If K is symmetric positive, then 

(v,Ku) + (v,Mu)g = (K*v,u) + (M*v,u)g (1.12) 

Proof - The proof follows from Green's Theorem. Ly definition we 


have 
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(v,Ku) - (K*v,u) = f v • (a 

a 


• Vu) + v • (V • (au) ) + v • Gu 



(v,Ku) - (K*v,u) =2 f n • (v • au) =2 f v • pu 
by Green's Theorem, and since p = n • a* Finally 

(v,Ku) - (K*v,u) = f (p*v •u + pv'u-v’M.u + v 1 0u) 

J bn 

= (M*v,u) b - (v,Mu) B 

which proves the lemma. 

The "first identity" can now be used to obtain Friedrichs 
"second identity." 

Lemma 1.2 If K is symmetric positive, then 

(u,Ku) + (u,Mu) b = (u,Gu) + (u,pu) B ( 1 . 13 ) 

Proof - It follows from the definitions of K* and M* that 
K + K* = G + G* and M + M* = p + p*; hence, letting v = u in 
the "first identity," we obtain 
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(u,Ku) + (u,Mu) B = i[(u, (K + K*)u) + (u, (M + M*)u)jQ 



G + G* 


2 



p + p* 
2 



= (u,Gu) + (u,pu) B 

The "second identity" immediately yields an inequality which 
will give us an a priori bound and insure uniqueness of any 
classical solution to a symmetric positive equation with semi- 
admissible boundary conditions. 

Lemma 1.3 Suppose u is a solution to (1.5) where M is 
semi -admissible. Let Aq be the smallest eigenvalue of 
(G + G*)/2 in iT. Then 


|u| 



(1.14) 


Proof - Since K is symmetric positive, Aq > 0, and therefore 
||u|| 2 S (u,Gu)/Aq. Using Lemma 1.2, since p + p* is non-negative 
definite by the assumption of the semi -admissible boundary condi- 
tion, we have 


ll u l | 2 s Ru,Gu) + (u,pu) B “l = -i. (u,Ku) , 

A G a g 

since Mu = 0, so that 

Hull 2 S i ||u|| ||Ku|| = i. ||u|| llfll 

One other inequality can be obtained for ||u|| B by as sinning 
that p + p* is positive definite. 

Lemma 1.4 Let u satisfy equation (1.5) where M is semi- 
admissible. Further, assume that (p + p*)/2 is positive definite 
on Sn with smallest eigenvalue A . 

H* 


Then 
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IMIb = —=■ Ilf II (1-15) 

VAgV 

Proof - From the hypothesis, 

IMlf S _L (u,|iu) S A [(u,nu) + (u,Gu)~] = A (u,Ku) 

^ V V V 

< f Hull ||Ku|| < _i_ ||f|| 2 

V V'G 

by Lemma 1.3. 

1.3 Uniqueness of a Cj_ Solution 

Lemma 1.3 insures the uniqueness of a classical solution, and 

2 

also that it is well posed in L for homogeneous boundary condi- 
tions . 

Theorem 1.1 If ueC-^fi) satisfies equation (1.5) where M is 
semi -admissible, then u is the unique solution to (1.5). Further 
(1.5) is well posed in the sense that for any e > 0 there exists 
a 6 > 0 such that if f is replaced by f e in (1.5) with 
|| f e - f || < 5, and if a solution Ug still exists, then Hug - u|| < e. 
Proof - Suppose that veC^fl) is any solution of (1.5), then 
K(u - v) = 0, M(u - v) = 0 is semi -admissible and by Lemma 1.3, 

||u - v|| = 0. For the second part let 6 = A q€, then 
K(ug - u) = fg - f, M(u€ - u) = 0, 

hence 

||ue - u|| £ A ||f 6 - f|| < e 
A G 

Actually piecewise differentiability of u is adequate for 
the above theorem provided u is continuous. This follows easily 
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since, when Greenes theorem is applied, the values of u along the 
discontinuities of the derivative will cancel, providing us with 
all the previous results . 

1.4 Weak and Strong Solutions 

By widening the class of solutions to (1.5) to include weak 
solutions it is quite easy to prove existence of a solution to a 
symmetric positive equation under only semi -admissible boundary 
conditions. We will use Friedrichs' definition of weak solution. 

Let V = Cp(ft) P|{v|M*v = 0 on A function ueH (defined in 

section 1.2) is a weak solution of (l.5) if fsH and for all v^V 

(v,f) = (K*v,u) (1.16) 

It follows from the "first identity" (1.12) that a classical solu- 
tion is also a weak solution. 

Theorem 1.2 If M is semi -admissible, there exists a weak solution 
to (1.5). 

Proof - Let -?^be the subspace of all functions w, where w = K*v 
with v€ V. Since K* is symmetric positive and M* is semi- 
admissible, Theorem 1.1 implies that v is unique for any given 
w. Hence, for any fixed feH, we can define a linear functional 
Lf, defined on<Pf^C by 

Lf(w) = (v,f). 

This linear functional is bounded, since 

|(v,f)| S Ml ||f|| < -L||f|| INI 

*o 

by Lemma 1.3 applied to K* and M*. By the Hahn -Banach theorem 
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Lf can be extended to all of <H, and by the Riesz representation 
theorem there is a ueH such that 

(v,f) = (wjii) 

which proves the theorem. 

This only shows that ueH, however, if ueCp(ft), we see from 
Lemma 1.1 that 

(v,Ku) + (v,Mu)p = (K*v,u) + (M*v,u) 

= (v,f) for all v€V. 

Hence (v,Ku - f) = 0 if v = 0 on dft, so that we must have 
Ku = f in ft. This in turn shows that (v,Mu)g must be zero. 
Friedrichs [l] shows that if, in addition, M is admissible, then 
Mu = 0. The conclusion then is that a weak solution which satis- 
fied admissible boundary conditions and is continuously differentia- 
able is also a classical solution to (1.5). 

A function ueM is a strong solution to (l.5) if there exists 
a sequence V) of functions such that each u^eCp(ft) and 

i3 ui - u ii + ii Kul - f ii + imi^} = o 

Variations of the definitions of weak and strong solutions are 
common (cf. Sarason [3]). In general it is not known whether a 
weak solution is differentiable; it is, however, possible, under 
certain additional hypotheses, to show that a weak solution is also 
a strong solution. One hypothesis used by Friedrichs [l] is that 
btt has a continuous normal. Sarason [3] considers the case where 
bQ is of class Cg- Sarason also considers the two-dimensional 
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case with corners, which requires special conditions to he satis- 
fied at the corners. Other "weak= strong" theorems are given in 
Sarason [3], Lax and Phillips [4], and Phillips and Sarason [5]. 

1.5 A Simple Example 

An illustration of the types of boundary conditions with more 
or less boundary data than usual can be given by means of a one- 
dimensional example. Suppose that 

Ku = 2x — + 2u = 0 for -1 S x £ 1 (1.17) 

dx 

If we write K in self adjoint form 

Ku = x — + + u 

dx dx 

we have a = x and 6 = 1, so that K is positive symmetric. At 
x = -1, 3 = na = -x, and we can let p = | 3 | = -x. Hence 
M = p - 3 = 0 and no boundary condition is imposed at x = -1. 

At x = 1, 3 = x, and letting p = |3|, we have again that M = 0, 
and no boundary condition is necessary at the right end either. 

Thus, for equation (1.17), no boundary condition at all is an 
admissible boundary condition 1 . To see that this is so, we can 
calculate the solution to (1.17). Since Ku = 2 d(xu)/dx = 0, we 
have xu = c, as the general solution. However, the theory is con- 
cerned only with solutions in L^(-l,l), and u = c/x is square 
integrable only for c = 0, so we do indeed have a unique solution 

r> 

in L (-1,1) without specifying any boundary data at all. 

A simple example can also be given of an ordinary differential 
equation which requires more boundary data than usual. For this let 
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Ku 


-2x 


du 

dx 


-2 


( 1 . 18 ) 


In self adjoint form 


Ku = -x fi - + u 

dx dx 

so that a = -x and G = 1. In this case if we make p = |p|, 
we get p = -0, so that M = p - 0 = 2 , at both x = 1, and 
x = -1. Hence, boundary data must be specified at both end 
points for admissible boundary conditions. Again, we can check 
this by solving the equation. The general solution to (1.18) is 

u = log |x| + c 


Since 



we see that we have a valid solution for 


any c. Also, because of the singularity at x = 0, we can 
specify the value of u at both x = 1 and x = -1. 


CHAPTER II 


FINITE DIFFERENCE SOLUTION OF SYMMETRIC POSITIVE 
DIFFERENTIAL EQUATIONS 

2 .1 Finite Difference Approximation to K and M 

First we will express K in a form slightly different from 
(1.3), by the use of (l.2). We have 

Ku = a • Vu + V • (au) + Gu 

= 2 V • (au) - (V • a) u + Gu (2.l) 

Using the concept of vectors whose components are themselves 
matrices or vectors leads to somewhat simpler notation for the 
application of Green's theorem. 

Lemma 2.1 (Green's Theorem) Let g be a continuously differentia- 
ble m-dimensional vector -valued function defined on A C R m , with 
vector components in either R, R r or R r X R r . Then 



Proof - Consider the case when g has matrix components, i.e., 
g =. (g'Sg 2 , • • ■> g m ) where g 1 = (g^ k ) is an r X r matrix. 


Then 
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is a matrix. Using the subscript j,k to indicate the element in 
the row and k^* 1 column, we have 



(using obvious notation ) ; therefore 



Similarly, the result holds when g has vector components, so 
the lemma is proved. 

We now integrate the equation Ku ** f over any region P Q !) 
using (2.1) and Green's theorem to obtain 




= 2 f au • n - f (V * a)u + J Gu *■* 
"dP J F J V 


= 2/* 0u - f (7 • a)u + f Gu = f f 
J bP «T J P J P 


(2.3) 


By a suitable approximation to (2.3) the desired finite difference 
equations will be obtained. 

Let H be a set of N mesh points for ft. It is not required 
for the theory that the mesh points all lie in ft. With each mesh 
point we identify a mesh region, Pj C ft by 



x kl * Vx k eH > k J* Ji 
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If Pj is adjacent to we say that Xj is connected to x^ 

(corresponding to the fact that the directed graph of the resulting 
matrix will have a directed path in both directions between j and 


k. 

see p. 16, [ 6] ) . 

Let 

Z j,k = l x j ' x kl' 

where 

x j 

is connected 

to 

x^, and let h 

= max 

^j,k' ^ow & e; fi ne 

A i 

to 

be 

the "volume" 

of 

Pj 811(1 L j,k 

to be 

the "area" of the 

r - 

1 

dimensional 


surface" between Pj and P^. We put Tj ^ = Pj P) P^* Figure 1 


illustrates mesh points and corresponding mesh regions for two 
dimensions. This concept of mesh regions is based on the sugges- 
tions of MacNeal [7]. We will always use the notation ^ to indi- 

j 

cate a sum over all points, x,, in H, and to indicate a sum 

over points, x k , which are connected to some one point, xj. 

The desired finite difference equation can now be obtained by 
a suitable approximation to equation (2.3). We use the symbol == 
to indicate the discrete approximation that will be used for each 



where Uj = u(x^) and p. v is the value of p for P.* at the 
center of rj^. (Note that Pj^k = “ Pk j)' ® ie approxiroation 
to the next term of equation (2.3) requires approximating u with 


Uj first, and then applying Green’s theorem before approximating 
a. With this we obtain 
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The final approximation is then 


(2.5) 


L Pu J * L J,k3j,k u J (2 ' 6 > 

r j,k 

Equations (2.4) and (2.6) take care of the integration over the 
interface between any Pj and P^. Now we need to make an approxi- 
mation for the boundary sides. It will be convenient to be able 
to subdivide P., Pi cJft into more than one piece. We will label 

each piece r j,B and we will use the convention that will 

B 

mean a summation over the B for just one j . We use l . - to 

3 

denote the distance from xj to xg, where xg is located at the 
" center" of r kj>B use< ^ B° r Bhe area of T^g. 

Also Pj^b = P( X B^ * Bhis notation is indicated for the two- 
dimensional case in Figure 1. The desired approximations are now 
given by 


L PU B L j , B^ j , B U B 


j,B 


f P u j T L j,BPj,B u j 


(2.7) 


( 2 . 8 ) 


j,B 


Finally the remaining terms in equation (2.3) are approximated by 


r a *- k ) a ) “j 


,p j 


(2.9) 
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L 


f = A i f 1 
• J J 


( 2 . 10 ) 


where Gj = G(xj ) and f j = f(xj). Also we can approximate J“ Ku 
by 

f Ku=A <j (K h u) j (2.11) 

P j 

where is the finite difference operator to be defined and 

which will approximate K. Using approximations (2.4) to (2.11) 
in equation (2.3) we arrive at the following definition of 

V K h u) j = X + u k) + 2 2 L j 

k B 

- 2 L j,kPj,k u J ' S + A j G j u j 


L j,kPj,k u k + 2 L j,B3j,B( 2u B - u j) + A J G J U , 

B 


j (2.12) 


where u here denotes a discrete function defined on H = H 
and uj = u(xj). We will seek to find a function defined on H 
and satisfying (K^u) j = fj for every x^eH. Of course the solution 
is not yet uniquely determined since there are more unknowns than 
equations. The boundary condition Mu = 0 will furnish us with 
the necessary information to determine u uniquely on H (but not 
necessarily on all of H) . 

Using to denote the boundary operator used to approximate 

M, we make the following definition 

(^)j,B = - J3j, B ( 2u B - V 


(2.13) 
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for all j where Pj is a boundary polygon, and for all boundary 
surfaces of P* (each of which is associated with a point x*) • It 
is easily seen that M h is consistent with M (i.e., (M^u) j,B 
Mu(xj g) as h -* 0 if u is continuous). The reason for this 

choice of is that the condition M^u = 0 can be used to 

el imina te ug in K^u in a simple manner, and also we will be able 

to prove basic identities for the finite difference operators 
analogous to those for the continuous operators (eqs. (1.12) 
and (1.13)). 

2.2 Basic Identities for the Finite Difference Operators 

The existence and uniqueness of a solution to the finite 
difference equation and the convergence to a continuous solution 
as h -*■ 0 depends on proving the basic identities for the dis- 
crete operators. Let be the finite dimensional Hilbert space 

of discrete functions defined on H. The inner product is given by 


(u,v) h = 



Vj,Xj eH 


(2.14) 


and 


Il u ll h = (u,u) h 

Also a "boundary" inner product is given by 

I 

j B 


(2.15) 


- 2 2 


B ' V j,B 


(2.16) 


for P.: a boundary mesh region, and 

J 

Hull? = (u,u) 


B. 


(2.17) 


22 


discrete adjoint operators K* and M* are defiined in 


the obvious way. 




(M h u) j,B - U J,B U J +B J,B (2U B - U j> (2 - 19) 

We can now give the "first identity" for the discrete operators. 
Lemma 2.2 If K is symmetric positive, then 

(v ,^h u) h + " (K h^ u) h + ( ^ u) Bh (2 ’ 20) 

for any functions u,v defined on H. 

Proof - Using the definitions, equations (2.12) and (2.18), we have 




By rearrangement, since 3j,k = “Pj,k> an< ^ s i nce Pj,k i s symmetric 


we have 
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2 



L j , kP j , k v k 



u k= -ZZ L j,k V j‘ ^ j , k u k 

J k 


and we see that all terms cancel with the exception of the boundary 
terms, so that 

(T »V>h - (I ^’ u) h "22 ' ^,b (2u b - u j> 

J -B v 

+ <\j,B (2V B - T j> • U j) <2 - 21 > 
On the other hand, using equations (2.13) and (2.19) 

- (v.VV' 2 2 L J,b(^,B V o ' U 0 + 2 J,B< 2l 'B-’'j>- »j) 

- ^ 2 Lj ' B f 3 ’ ' T j ‘ p J.b (2u b - u j>) 

•which is the same as the right side of (2.21). Hence the "first 
identity" for the difference operators is proved. 

+ K £ = 

G + G* and + M* = \i + p*. By letting v = u in (2.20) we 
can prove the discrete " second identity 11 exactly as for the con- 
tinuous case (Lemma 1.2). 

Lemma 2.5 If K is symmetric positive, then 

(u,K h u) h + (u^u)^ = (u,Gu) h + ( u, |iu) (2.22) 

2.5 Existence of Solution to Finite Difference Equations 

Using equation (2.13) and M^u = 0 we can eliminate Ug from 
equation (2.12) so that the equation K^u = f can be reduced to 


The discrete operators have been defined so that 


24 


2 L «j,kPj,k u k + ^ L j,B^j,B u j + A j G j u j “ A j f j^ V j (2.23) 
k B 

If we consider the case when Cl is two-dimensional and rectangular, 

and the P. are all equal rectangles, we can compare (2.23) with 
J 

the finite difference equation obtained by Chu [2]. The equation 

obtained by Chu is the same as (2.23) for interior rectangles, but 

is different for boundary rectangles. 

Let A be the rN X rN matrix of coefficients of (2.23). 

Letting (u,v) = ,£ u * v , the ordinary vector inner product, we 

j 0 j 

have 

<u,Au) = (u,K h u) h + (u,!^)* (2.24) 

Hence, by the ’’second identity” (2.22), A has positive 
definite symmetric part which shows that A is non -singular. We 
can also obtain an a priori bound for ||u||^ just as in the con- 
tinuous case. 

Lemma 2.4 Suppose u is a solution to 

K^u =5 f, M^u = 0 

where K is symmetric positive and M is semi -admissible. Then 

Nlh S ^7 ll f llh (2-25) 

If in addition, (p + p*) is positive definite on dn, then 


u 


VA g A u 


M h 


(2.26) 
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Proof - The proof is identical to that for Lemmas 1.3- and 1.4, 
but using the h norms and inner products. 

2.4 Convergence of the Finite Difference Solution to a 
Continuous Solution 

It is possible to show that the solution of the finite differ- 
ence equation (2.23) converges strongly to a continuously differ- 
entiable solution of equation (1.5), under the proper hypotheses. 

For simplicity we prove convergence only for the case when 0 is 
two-dimensional (m = 2). Extension to regions in higher dimen- 
sions, with the same rate of convergence, follows directly. To 
allow the type of comparison we wish to make we will define 
operators mapping 4 into and vice versa. Let r^: £{ -* 

be "the projection defined by 

(r.u). = u(x.) for all x.€H (2.27) 

J J J 

In the other direction, let p^: -*• be an injection mapping 

defined by 

p h u h^ = ( u h^ j > for a11 xeP j (2.28) 

We immediately have the following relations, 

r^p^ = I (2.29) 

llp h u hH = IMh for a11 u lA (2.30) 

We can now state our basic convergence theorem for two-dimensional 
regions . 

Theorem 2 . 1 Suppose that u€C^(0) satisfies 

Ku = f on a Q R 2 
Mu = 0 on So 
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where K is symmetric positive, and p + p* is positive definite 
on bSi. For any given h > 0, let H h be a set of associated 
mesh points such that the maximum distance between connected 
nodes is less than h and also that Lj^, anc * l x “ x jl 

for xePj are all less than h. It is assumed that the mesh is 
sufficiently regular so that h^/Aj for each Pj is bounded 
independently of h by a constant Kp > 0, which is possible for 
sufficiently nice regions. Also it is assumed that a uniform 
rectangular mesh is used for all Pj any point of which is at a 
distance greater than Kgh from c>Q, where Kg is a positive 
constant. It is assumed that ot£C^(fi). 

Let Uft^h be the unique solution to 

K h u h = r h f on H h 

Vh = 0 


Then Up^u^ “ u ll = 0(h V ) as h ■* 0 for any positive v < l/2.| 
Chu [2] proved convergence of his finite difference scheme, 
where Q is a rectangle or a region with four corners, but the 
rate of convergence was not established. 

Proof - Define w h = u h ” r h u ‘ ^et ^'* le sma ll es t eigen- 

value of (G ,+ ;G*)/2 in ft. Using the "second identity" (2.22), 
we have 


l|w h ||g£-L j^(w h ,Gw h ) h +(w h ,pw h ) Bh » i. [Iw h ,K h w h ) h +(w h ,M h w h ) Bh ] 


Using the Cauchy-Schwartz inequality, we have 
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Ml ^ ^ ( H w hllh HVhHh + \K\ HvJb^ (2 * 31) 

We will show that ||K h w h || h = C^h 1 / 2 ) and IN^hllB^ = as h 

We shall need the following lemma. 

Lemma 2.5 Let g be a function defined on a finite region PC^ 2 > 

and suppose that g satisfies a Lipschitz condition, i.e., there 

is a constant K3 > 0 such that |g(x) - g(y) | ^ %| x - y|, 

for all x, yeP. Then, if A 0 is the area of P and |x-x 0 | < h 
in P, 


S(x 0 ) - ^ J g(x) 


^K 3 h 


Proof - By direct calculation 


g(x Q )-i j (gt xB ~ T 2 f fs(x Q ) -g(x)) 




We proceed now with the proof of the theorem. Let 0^ denote 

that portion of Q consisting of those P* which are rectangular, 

j 

and let Qg denote the rest of the Pj. From the hypothesis we 
see that the area of (1% is less than the length of <5Q times 
K^h. We have now that 

ll K h w hlih= / (Ph K h w h ) 2 + f (Ph K h w h ) 2 

“h 2 

= / , f ( Ku ( x j ) -(K h r h u)j) 2 + 2> I ■ vlv h i h u ' , j 


p. 

J 


2 + ^ / (Ku(x.) - (K^r^u) J 2 

0 eJ 2 p 


(2.32) 


where 


J i = {JlPjCaJ, i = 1,2 
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To simplify notation we will use Uj for u(xj) and u-g 
We now obtain a suitable bound for |Ku(xj) - ( K-^r^u ) j | 


| Ku(xj ) - (K h r h u)j| = 1 2V * (au)(xj) - (V-a)u(xj) + Gju-j 


-2% 


Pj,k (u j + u k> - 2 


"a* 2 3 J,b u b 

B J 




B 


A” ®J,B“3 - °j u j 
J 


2V • ( a u) (xj ) - 2 ~Aj~ 0j,k<V U k )_2 Z ~Aj^ 
k B 


(V-a)u(xj)-2^e J)k u r 2^ Pj,B“j 
k B J 


Consider the first term in the last expression above 

2V- (au)(xj) - y -^Bj, k (uj+u k )-2 y ^ P 3 , B U B 
k 0 B ° 


2V • (au)(x.) -T- f V • (au) 

0 Aj y 

X / 2 ^ u -<p">j,k) 


j,k 



for u(x B ). 


p j,b u b| 

(2.33) 


s j,k < 2 u j,k- (u j +u k )) 


(2.34) 
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By Lemma 2.5, since a and ueC 2 (n) imply that their derivatives 
satisfy a Lipschitz condition, 


2 V 


(au)(xj) - Y; f V'(au) 


3 P J 


= 0(h) 


(2.35) 


We consider now the case when so that Pj is a 

rectangle with Xj at the center. 


Since uec 2 (ft), we have 


- ¥ 


ui X + l J£ u" 


j> k Ttjz 


<*L> 


^ = + "2^ U J,k + ¥ U " (lz) 

where the derivatives are directional derivatives in the direction 
x k " x j* Hence, if |u"| < K g in Q, we have 


K 


l 2u j,k - + u kM < i 


^h 2 


This means that 


f -< u j +u k 

Jt. , 


)) 


j jk 


^ k j,k IIPj,kll l 2u j,k “ + u k^ l = ) 

(2.36) 


when je J-j-. 

We now examine a Taylor series expansion for pu about the 
point x^ k =(xj + x k )/2. 
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' 2 ~ 
3(x ja + tz)u(x J#k + tz) = (pu)j^ + t(^ (pu)) J#k + \ g(|l) 

2 (2 

* (x J,k - tz)u(x Jjk - tz) = (pu) j>k - t(A (pu)) J#k + V g(|2) ^ 

where z is a unit vector orthogonal to Xj - x k , t is a scalar 
parameter, g(|) = (g 1 (l 1 ),g 2 (i 2 )^ • ■ •> g r (l r )), g± is the i th 
component of the vector d 2 /dt^ (pu), and is a point on the 

straight line between xj^ k + (L^ k /2)z and xj^ k - (Lj^ k /2)z t 
Using (2.37) we obtain the following bound, 


f eu-Ou) j(] 


L i .k / 2 




(3(xj^ k + tz)u(xj k +tz) 


+ P^ x j,k " tz ^ u ^ x j,k ' tz ) _ 2 (^ u )j >k ) <3 -t 

£ f t 2 K4 dt = 0(h 3 ) 


Now, using (2.35), (2.36) and (2.38) in (2.34) we obtain 


|2V*(au)(x.) - V 3j k (u + u k )| = 0(h) (2.39) 

k J ' 

for all jsJj, since h 2 /Aj and the boundary terms are not 

present. 

Consider now the second term on 'the right of (2.33): 
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(V»a)u(* 3 ) - ^ ^ - X S J,B“ 

k J B 3 

1 


(V • a)u(xj) - / (V-a)u + '£7 f (V*a)(u~Uj) 

.It-* J 


J P 0 


P J 


y / (p-3j,k) u j + y / (p-p^b)^ 

* ' ^ Tl -p _ 

B P j,B 


k P j,k 


(2.40) 


By Lemma 2,5 


(V • a)u(xj) - J- J (V * a)u 


p o 


= 0(h) 


(2.41) 


Next, since u satisfies a Lipschitz condition, |x - x j | < h for 


all xePy, and since ||v • a|| is uniformly bounded in Q, we have 
J 


a j r p j 


(V • a)(u - u^ 


= 0(h) 


(2.42) 


Finally, since p. , and p. ^ are each evaluated at the midpoint 

Jj-U 

of or respectively, we can use a Taylor series 

analysis, as in deriving equation (2.38), to obtain 


JL_ 

A. 


Y / 

J p 


k 


(3 " 3j,k) u j + y f_ (3-3j,B)^j 




J T 

B ^ <5 >B 


= 0(h) (2.43) 


Combining (2.41), (2.42), and (2.43) in (2.40) we obtain 


(V-a)u( Xj )-2^ 3j, k Uj " 2 ^ 3j,B^j 
k J B J 


= 0(h) 


(2.44) 
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Note that (2.44) holds for all j, not just for jeJp . 

We can now substitute (2.39) and (2.44) in (2.33) to obtain 

|Ku(xj) - (K h r h u)j| = 0(h) for all jej^ (2.45) 

We cannot obtain as good a bound for |Ku(x^) - (Kh r h u )jl 
when jejg, although (2.44) holds, since ^ is not in general 
bisected by the line between xj and x^. However, we can show 
that |Ku(xj) - (K^r^u) j | is uniformly bounded for je J 2 , which 
is adequate since the area of U 2 is of order h. The two in- 
equalities which must be re-examined are (2.36) and (2.38). 

We now have, since u and ((3u) satisfy Lipschitz conditions, that 


/ Pj,k( 2 u j,k'( u j + u k)) 


j,k 


= 0(h 2 ) 


(2.46) 


/ 

= 0(h 2 ) ^ 

r j,k 


f - (3 u )j,B 

= 0(h 2 ) 

r j,B 

J 


(2.47) 


Using this, with the other results which still hold, we see that 
|Ku(xj) - (K h r h u)j| is uniformly bounded for jeJ 2 , as h -*• 0. 
Using this, together with (2.45) in (2.32) we obtain 


H K h w hHh = o(h ) + 0(h) 


(2.48) 


so that 


||K h w h || h - CKh 1 / 2 ) 

The next step is to show that = 0(h). We have 

II^JIb = H M h u h " M h r h u llB h = H M h r h U llB h 


-:>■ ' 


(2.49) 
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since M^u^ = 0. Now 

ll M h r h u ll?h = X 2 L j,B( M h r h u )j,B 

J B 

" 2 X L j,B(^j,B u j - 0 J,b( 2u B - u j)) 2 

However, using the fact that Pj } b u b = 1_l j > B u B^ 

h,B u j -3j,B( 2u B“ u j)| ^ l^j,B( u j- u B)| + |Pj,B(^B- u j)| = °( h ) 
since u is differentiable, and ||p|| and ||p|| are uniformly- 
bounded. This shows that 

iiwnib = °<» 2 >> 


since L. t> is simply the length of cto. This proves that 

ifs J,B 

\K*h\ = 0(h) 

Using (2.49) and (2.50) in (2.3l), we see that 

Kllg = l|v h || h 0(hl/2) + KHB h 0 ( h > 

From Lemma 2 .4, ||wjJ|g^ must be bounded, since 

+ H r h u llB h 


(2.50) 


(2.51) 


■ vrr l|rhf|lh + l|rhU|lB h 

^ Gp 

which is certainly uniformly bounded as h -*• 0. Likewise ||wjJ| h is 
bounded. So from (2. 51) we have 

||v h || h = 0(h!/ 4 ) (2.52) 

However, if we use (2.52) in (2.51) we get llv^llft = 0(h^/®), or by 
repeating this procedure enough times, 

||w h || h = 0(h v ), for any positive v < l/2 


(2.53) 
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Finally, we establish the convergence rate for - u|| . 

We have 

l|Ph u h ' U H ^ HPh u h ~ Ph r h u ll + HPh r h u ' U H 

= l|w h || h + ||p h r h u - u|| (2.54) 

The last term can be estimated by 

l|p h r h u " U H 2 * 2 f " u ^ 2 = °^ h2 ^ (2.55) 

J p j 

Using (2.53) and (2.55) in (2.54) we get 

||p h u h -ull =0(h v ) + 0(h) = 0(h v ), for any positive v < l/2 (2.56) 

This completes the proof of Theorem (2.1). 

2.5 Solution of the Finite Difference Equation 

For our method to be of practical use we must have some 
method for computing the solution to the finite difference equa- 
tion (2.23). We will consider only the two-dimensional case 
here. In any case we can partition the matrix A so as to be block 
tridiagonal. For example, suppose that the mesh points H lie on 
lines such that the mesh points on any one line are connected 
only to points on the same line or adjacent lines. Then we can 
partition A into blocks corresponding to each line. The diagonal 
blocks will themselves be block tridiagonal with r X r blocks. 

The matrix equation can then be solved by the block tridiagonal 
algorithm ([8] and [6], p. 196). We suppose A to be written in 
the form. 
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where NL is the number of lines. Each B.^ is an rn X rn 
block tridiagonal matrix, where n is the number of points on 

-Ll_ 

the i n line. From equation (2.23) since 3^ k = -£ k j we see 
that = C^_^. Thus need not be stored for a computer 

solution. The block tridiagonal algorithm is completely analogous 
to the ordinary tridiagonal algorithm. Suppose the equation to be 
solved is Au = f, where u and f are partitioned, as required. 

A typical block equation is 

^i u i-l + ®i u i + G i u i+1 = ^i 

First let 


w i - 


>i- f i 


The forward sweep is given by 


G. = a.w: 1 .. 

l li-l 


y i ~ f i " Vi -3 


W i = B i - G i C i-l 


/ for i = 2,3, 


J 


NL 


This is followed by the backward sweep. First, 
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U NL = W fji y KL 

Then 

u i = ^(yi " C i u i+l) for i = NL - 1, NL - 2, . . . , 1 

Of course this algorithm will not work for every non-singular 
block tridiagonal matrix. However, Schecter [8], gives a suffi- 
cient condition for the validity of the algorithm, and that is 
simply that A has definite symmetric part. We have already shown 
that A has positive definite symmetric part. There is one real 
disadvantage to the method, however, and that is the fact that 
each WJ^" is a full matrix and must be stored during the forward 
sweep for use on the backward sweep. This results in large com- 
puter storage requirements, and the use of tapes or disks for 
temporary storage for only a moderate number of mesh points. This, 
of course, is very time consuming. An alternate procedure is 
suggested by Schecter [8]. In Schecter' s method only one matrix 
need be inverted and stored for a number of consecutive lines 
with an equal number of points per line. However, the matrix 
to be inverted may be ill-conditioned if too many lines are grouped 
in this way. 

An alternate method of solution may be possible in some cases. 
Note that A may be decomposed as 

A = D + S 

where D is Hermitian and positive definite, and S is skew 
symmetric. The eigenvalues of D are usually easy to calculate 
since D is block diagonal with r x r blocks. If the smallest 
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eigenvalue, Ap, of D is larger than the spectral radius, p(S), 
of S, we will have 

!|d - 1 s|| ^ IlD- 1 !! |}sj| = i p(s) < 1 

In this case we could use the following iterative method. Let 
u^ 0 ) be arbitrary, and define u^ recursively by 

Du^ 1 ) = -Su^ 1 ' 1 ) + f 

In this case lim u^ = u. In general, though, the eigenvalues of 

i-wo 

D will not all be sufficiently large for this simple method to 
work. However, the original finite difference equations can be 
modified in some cases by the addition of a "viscosity" term, so 
as obtain a convergent iterative procedure for the solution of the 
matrix equation. This will be discussed further in Chapter III. 

2.6 Convergence to a Weak Solution 

We can consider the discrete analogue of a weak solution. Let 
be the set of discrete functions, v^, defined on H and 
satisfying M^v^ = 0. For a discrete weak solution, u^, we would 
then require that 

(Kjv h ,u h )h = (vh^ r h f )h for a11 V ^ v h (2.58) 

Form the "first identity" (2.20) we have then 

(v h> r h f) h = K'ViA + ( v h^ u h\ for a11 ^h (2 - 59) 

We see from this that (K^u^) ^ = f ^ for all Pj which are not on 
the boundary, by choosing (v^)j = 1, and (v^^. = 0 for k ^ j. 
Because of the discrete nature of the equations we are not assured 
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of satisfying the boundary conditions. However, conversely, 

if satisfies = r^f and M^u^ = 0 we see imeediately 

that (2.58) must be satisfied. 

Chu [2] has shown weak convergence of his finite difference 
solution to a weak solution of a symmetric positive equation and 
Cea [9] has investigated generally the question of weak or strong 
convergence of approximate solutions to weak solutions of elliptic 
equations. Using these ideas, we can prove weak convergence of our 
finite difference solutions to weak solutions of symmetric 
positive equations. 

Theorem 2.2 For any h > 0, let be a set of mesh points 

satisfying the requirements of Theorem 2.1. It is assumed that 

o _ 

q,eCr(f2). Let u^ be the unique solution to 

^h u h = r h^ 


Vh= 0 



is a positive sequence converging to zero, then 
has a subsequence which converges weakly in H 


to a 


weak solution, u, of equation (1.5), that is 

(K*v,u) = (v,f) for all veV 
Furthermore, if u is a unique weak solution. 


then 



converges weakly to u. 
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Prbof - First we note that ||Ph u hll is "bounded, since 

l|Ph u hll = ll u hllh ^ T~ il r h f Hh> Lemma 2.4. Hence, there is a sub- 

^AJ 

sequence of {Pfc . u h • converges weakly to some usU. (See 

Theorem 4.41-B, Taylor [10].) For convenience of notation we will 
suppress the subscripts on the h. 

We have, for all veV, 

|(Kjr h v,u h )h - (K*v,p h uj 1 ) | £ |Ph K h r h v " K * V H IIPh u hll 

£ (\K r h v ~ r h^ y Hh + HPh r h^ v " ^ v ll)llPh u hH (2.60) 

But Hp^r^K^v - K*v|| 0, and in Theorem 2.1 we can substitute K* 

for E in equation (2.49) to show that ~ r^K^vH ■* 0 

(since K^w^ = r^Ku - K^r^u) . Since IIp^u^H is bounded, 

lim | (K^r h v,u h ) h - (K^pj^) | = 0 
h->0 

However, since K*ve <H, we know that lim (K^Vjp^u^) = (K*v,u) 

h-»0 

We have shown, then, that 


lim (K* r h v,u h ) h = (K*v,u) , for all veV. 
h-+0 ~ l 


(2.61) 

The discrete "first identity", equation (2.20), gives 

(Kjr h v,Uh) h + (M*r h v,u h ) Bh = (r h v,K h u h ) h + (r^M^)^* (r h v,r h f) h 

(2.62) 

Hence 

I ( K h r h y , u h)h ~ ( r h v > r h f )hl ^ H^ r h v llB h H u hllB h (2.63) 

By Lemma 2.4 ||u h ||g^ <; ll r h f llh'/V^G\t which is bounded. 
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Also, the proof of equation (2.50) 
for all v€V, so that 


shows that lim 
h-0 


IlMVi 


®h 


lim | (Kjr h v,u h ) h - (r h v,r h f) h | = 0 
h-*0 


Further, it is obvious that 

lim (r h v,r h f) h = (v,f) 
h-*0 


Combining (2.61), (2.64) and (2.65) gives 

(K*v,u) = (v,f), for all v£V, 
which completes the proof of the theorem. 


= 0 , 

(2.64) 

(2.65) 



CHAPTER III 


SPECIAL FINITE DIFFERENCE SCHEME FOR ITERATIVE 
SOLUTION OF MATRIX EQUATION 

5.1 Special Finite Difference Scheme 

As pointed out in section 2.5, the matrix equation Au = f 
can be solved by an iterative procedure if the eigenvalues of the 
diagonal coefficient matrix are sufficiently large compared to the 
eigenvalues of the off-diagonal coefficient matrix. Following the 
idea of Chu [2], we modify the finite difference equation by adding 
a "viscosity" term which will have a diminishing effect on the fi- 
nite difference equations as h->-0, and yet will assure the conver- 
gence of an iterative method. Unfortunately, the method is not 
applicable to every arrangement of mesh points. In fact there are 
rather severe restrictions which must be met. The first require- 
ment is that the difference in areas of adjacent mesh regions be 
sufficiently small. This cannot be readily done along an irregular 
boundary, however, unless the boundary is modified. A problem 
arises if the boundary is modified. The boundary condition is 
given by Mu = (p - p)u = 0 on dfi. We need to extend M to be 
defined in a neighborhood of the boundary. It is possible to 
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extend M continuously in a neighborhood of the boundary. How- 
ever, if the direction of the boundary changes, 0 changes 
drastically, and we have no assurance that p will be positive 
definite. The second requirement then is that M can be extended 
continuously over a neighborhood of the boundary, in such a way 
that p will have positive definite symmetric part along the 
approximating boundary. 

Let £2^ be an approximation to £2. £2^ will have to meet 

several requirements to be specified later. H h will denote a set 
of mesh points associated with £2^ and with maximum distance h 
between connected nodes, and will denote U ®" ie 

discrete inner product is given by 

< u h> v h) = 2 A j (u h } j * (v h)j ( 3 -D 

3 


with the Aj being the area of Pj C £2^- Simiarly, the "boundary" 
inner product is changed so that the lengths, Lj^-g, are the lengths 
along d£2 h . 


We define now two new finite difference operators, and 
by 


(K, 


k B 


(3.2) 


oAj 


(M h U )j,B - (U J ‘ U B } 


(3.3) 
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where a is a positive number which must satisfy requirements to 
he •specified later. 

It will be useful to prove a slightly different version of the 
"second identity". 

Lemma 3.1 If K is symmetric positive, then 


2 

(u h ,K h u h ) h + (u h ,M h u h ) Bh = (u h ,Gu h ) h + (u h ,qu h ) Bh + —J- (uj -%)* 

' _ i 

)>k) 

(3.4) 


where / i indicates a sum over every (,i,k) pair where x^ is 

(j,k) 0 

connected to x k . 

Proof : Using the "second identity" for and M^, equation 

(2.22), we have 

(u h ,K h u h )h + ( u h , ^i u h)B h = ^ u h^ Gu h^h + ^h^h^ 



which completes the proof. 
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Lemma 3.1 immediately assures the existence and uniqueness of 
a solution for the special finite difference scheme. Using 
M^u^ = 0 to eliminate u-g from K^u^ = r^f , we obtain 



(3.5) 


for all x^eH^ 

Let A be the matrix of coefficients of (3.5). 

Lemma 3.2 If K is symmetric positive, then 

Vh * r h f 
= 0 

has a unique solution on H^. 

Proof : The hypothesis implies that 

(u,Au) = (u h ,K h u h ) h + (u^u^ 

By Lemma 3.1 A has positive definite symmetric part, and hence 
is non-singular. Thus (3.5) defines u^ uniquely on H^. 

Also it will be noted that the "second identity" of Lemma 3.1 
will give the same a priori bounds for Hu^H^ and ||u^||g^ as given 
by (2.25) and (2.26). 

5.2 Convergence of Special Finite Difference Scheme 

We will now show that the special finite difference scheme 
converges to a smooth solution, under a number of hypotheses 
given in the theorem. The theorem also includes all the hypotheses 
needed to assure convergence of the iterative matrix solution. 
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Though quite a number of requirements are given, there are only- 
two essential restrictions, namely, that the areas A. must be 
nearly uniform, and that M can be specified on a modified 
boundary in such a way that p remains positive definite. 

Theorem 5.1 Suppose that ueC^(n) satisfies 

Ku = f on f2 
Mu = 0 on 

where K is symmetric positive. For any h > 0, let be an 

approximation to fi, and let be a corresponding set of mesh 

points with maximum distance h between connected nodes, and 

also with and |x - Xj| for xeP^ all less than h. 

It is assumed that the following hypotheses are satisfied: 

(i) There exists Kq > 0, independent of h, such that for 

every P. we have h^/A . < K, . 

J J - 1 

(ii) There exists Kg > 0, independent of h, such that all 
Pj with any point at a distance greater than Kgh from dfi are 
equal rectangles. 

(iii) There exists Kj > 0, independent of h, such that for 
all xedn^, the distance from x to is less than K^h. 

(iv) There exists > 0, such that M can be extended so 
as to satisfy a uniform Lipschitz condition at all points at a 
distance less than from dft. 

(v) is such that n = M + g has positive definite 
symmetric part on 
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(vi) Let W be the set of points that are a distance less 
than K4 from df2. Then a, G, and f are all extended to be 
defined on with cc,eC 2 (f}(JW) and G positive definite on 

nUw. 

(vii) There exists K5 > 0, independent of h, such that all 
points, x i , associated with a boundary polygon, P., are in the 
polygon, and at a sufficient distance, 2^ g, from any boundary 
node, xg, of Pj so that Aj £ 

(viii) Either fyjCZ G or else u can be extended so that 
ueC 2 (JT h ) . 

(ix) 0 > rjK-. p-Q + d, where d > 0 and = sup p(n*a(x)), 
T B B xeftUw 

where n is any unit vector and t) is the maximum number of 
nodes connected to any one node. 

(x) |Aj/A k -l]< d7\ G (hf) 2 /(r) 2 a 2 h),. for all connected nodes, 

Xj and x k , where is the smallest eigenvalue of G in and 

h' = min(2^ k ) . 

(xi) The length of is uniformly bounded. 

Let u^ be the unique solution to 

Vh = r h f 

M h u h = 0 

K- r^ul^OChy) as hr*0, for any positive v < l/2 


then 
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Proof : Letting w h = u h ~ r h u ^ and - usln S "the second identity, 

(3.4), we see that the inequality (2.31) is still valid for 
and 

'IKJIh ^ ^l w hllh ll^h w hllh + l^hll^ ll^hHB^ ^ 3 ’ 6 ^ 

We have 

Vh = ^h u h - V h u= r h f -K h r h u = r fa f - K h r h u +K h r h u h -- K h r H u, 
hence 

llVhllh * ll r h Ku - K h r h u llh + ll K h r h u - V h u|| h (3.7) 

In checking the proof of Theorem 2.1 we see that r h Ku - K^r h u 
is the same as K h w h (Theorem 2.1), hence the bound of (2.49) 
holds for this term.> 

||r h Ku - K^ul^ = 0(h 1 / 2 ) (3.8) 

For the other term we have 

"(^ - ^ = Z (3 - 9) 

j \k B / 

Let J-, denote the set of subscripts for those P. which are 

equal rectangles, and let Jg denote the rest of the subscripts. 

When jeJ^ we have only the term ^(u^ - u^/Zj ^ to consider. 

k 

Because of the rectangular arrangement of points we can use a 
Taylor series analysis to show that 


u j - u k 


= 0(h) 
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so that 



= 0(h2) 


(3.10) 

J j 

j€ Ji V k / 

On the other hand, when jeJg we c 9 - 11110 ^ do as well. However, we 
note that both (uj -u k )/2j^ k and (uj - u B )/Zj^ k are uniformly 
bounded since u has a bounded derivative. Also, by hypothesis 
(ii), / , A, = 0(h), so that 

¥4 



u i_: u b | = 0(h) 


(3.11) 

v 4 J,b J 

jeJg ^k B 

It is assumed, of course, that the number of nodes connected to 
any one node is bounded as h+0. 

Now, using (3.10) and (3.1l) in (3.9) we have 

II (% - K h )r h u|| h = 0(h 2 / 2 ) (3.12) 

Taking this together with (3.8) in (3.7) finally 

||K h w h || h = Oth 1 / 2 ) (3.13) 

It is necessary how to obtain a bound for IlM^w^Hg • Since 

h 

Vh = Vh - Vh u “ -Vh u > we have 

ll% w hllB h ^ H M h r h u llB h + IK^h - 


(3.14) 


H M ti r h U llB X 2 L j,B ( ^j,B ’ P j,B (2U B ' 
.1 B 


)r 


We have 
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We can establish a bound, since 


b,B - »j,B< 2u B - u j>l £ b,B< u j - U B>I 

+ lb,B ’ + l^j,B^ u j ' VI 

The first and last term on the right are of order h, since u is 
differentiable and ||n|| and ||p|| are bounded. Ey hypothesis (iv) 

M satisfies a Lipschitz condition, and so does u. Since the 
distance from xg to do is less than K3I1 by (iii) and 

Mu = 0 on So, we see that |(nj b " B^ u bI = ^(h) * Since, by 

(xi) L is uniformly bounded, we have 

j B 

IlMh^ull^ = 0(h) (3.15) 


Also 

II («h - W4 



- U B> 2 ’ bj (Tll) 

j B 



= 0(h 2 ) 



(3. 

,16) 

This shows that 








H Mw 'hllB h = °( h ) 


(3. 

.17) 

We check now 

to see 

that ||w h || h and ||w h | 

L are bounded, 

\ 

We 

have, using the a 

priori 

bound for ||u h || h . 



# W A s 

K» h ♦ 

■ 11 r h u|l h s jf- i |r n f|l h 

+ l|r h uil h 

(3. 

.18) 
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which must he bounded since f and u are. In the same manner, 
||w h |L must be bounded. Using this fact together with (3.13) 
and (3.17) in (3.6) we have 

Kl| h - 0(h 1 / 4 ) (3.19) 

Using now (3.19) in (3.6) we get HwjJljj = 0(h^/®) and by repeating 
the process as many times as needed we get 

llwjlh = °( hV )> for any positive v < l/2 : (3.20) 

5.5 Convergence of the Matrix Iterative Solution 

For the iterative solution of the matrix equation Au = f 
we will split A into a block diagonal part D, and off diagonal 
part B. (We will suppress the subscript h on the finite 
difference solution u^. ) Thus, from (3.5), the j’*'* 1 block of D 
is an r X r matrix. 



I + 





and a typical block element of B is 

a A.* 

B j>k = L j,k P j,k " l I 

<] >&■ 

and A = D + B. The iterative method is given by 

u^ i+1 ^ = -D^Bu^ + D -1 f 

where u^) is arbitrary. The hypotheses of Theorem 3.1 assure 
the convergence of u^ 4 ^ to Ui 
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Theorem 5.2 For any h > 0, let and satisfy the hypotheses 

of Theorem 3.1. Let u^ be an arbitrary vector defined on 

and let ( u(i) )i=0 be a sequence defined recursively by 

u ( i+1 ) = -D -1 Bu( 1 ) * D -1 f 

Then lim u^^ = u, where Au = f . 
i-*<» 

Proof - By the contraction mapping theorem it is sufficient to show 
that ||D -b B|| < 1 for some matrix norm. Let v be an arbitrary 
vector defined on H h , and let w = D -1 Bv. Since Dw = Bv, we 
have 

(w,Dw) = (w,Bv) 


or 



This last inequality follows from the fact that 

(w,Hv> < i (w,Hw) + | <v,Hv> 
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for any positive definite Hermitian matrix. We see that 
(aAj)/(2j^ k ) I - Lj v^j,k is positive definite, since 
aA. aA. 

I . ^ h > K7 h > ^ L j,kP(Pj,k) (3.22) 

J . K _L 


by (i) and (ix). By rearranging the terms of (3.2l) so as to have 
all the w terms on the left and all the v terms to the right, 
we obtain 



The last expression was obtained by interchanging j and k, 
since 

2 j,k = \,y L j,k = \y and P j,k = ' P k,j 


We can write (3.23) in the following form. 
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or 

(w,Xw) + (w,Yw) ^ (v,Yv) + (v,Zv) (3.25) 

vhere x, Y, and Z are matrices defined by (3.24). 

We have already shown that Y is positive definite (using 
(3.22)); hence we can define a norm by 

IMIy = <v,Yv> (3.26) 

We will show that D”^B is a strict contraction in the Y norm. 
First we will need some inequalities. We have 


<w,xv> > (3.27) 

Next, we have 



by (i) and (ix), so that 

(w,Yw) < 12 . ||w||^ (3.28) 


Also (w,Yw) can be bounded below, since 



by (i) and (ix) . • Thus, we have 


< v > Yv > ^ ^ IMI^ (3.29) 


Finally, since 



< 
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we have 


<v,Zv> £ A M || v || h (3.30) 

where A = max j A^./A ^ - l| , for all connected nodes, x^ and x^. 
From the definition (3.26), and using (3.27) and (3.28) we have 

Acjh’ 


(w,Xw) + (w,Yw) > ^1 + 
On the other hand from (3.29) and (3.30) 
( v,Yv) + (v,Zv) £ 


Tier 


w 


i (£)~ 


(3.31) 


(3.32) 


Substituting (3,31) and (3.32) in (3.25) we have 


Mly < 


1 + !1^ /'JlY 


1 + 


r\a 



(3.33) 


Since w = D~^Bv, and v is arbitrary, we see that ||D _ ^B||y < 1 
since 

<%h' / h , 


A < 


2 2 
T) a 


(v) 


(3.34) 


by hypothesis (x) . This completes the proof of Theorem 3.2. 

Of course, if can be selected so that all the Aj are 

equal, then hypothesis (x) is satisfied, and 


I|d _ 1 b|Iy < 



In the special case where all the P. are equal rectangles, 
r) = 4, so that 


(3.35) 



v,.W 2 " 


( 3 . 36 ) 



CHAPTER IV 


APPLICATION TO THE TRICOMI EQUATION 

4.1 Transonic Gas Dynamics Problem 

An example of a problem of physical significance can be 
drawn from the field of gas dynamics. A stream function is 
introduced such that derivatives of the stream function are 
velocities of the gas. The stream function satisfies a second 
order partial differential equation which is elliptic where the 
flow is subsonic, and hyperbolic where the flow is supersonic. 

The equation is of mixed type, then, for a transonic flow problem. 
When the equation is linearized by means of a hodagraph trans- 
formation, and after a further transformation, the Tricomi equa- 
tion results, 

F(y)V - f =0 (4.1) 

xx yy 

where F(y) is a continuous monotone function such that yF(y) > 0 
for y ^ 0. Details of the derivation of (4.1) are given by 
Bers [ll] . A sqlution for (4.1) for a region Q which includes 
a portion of the x-axis is determined by proper boundary value 
data along portions of The proper boundary value data is 

known only for special cases. Usually the appropriate boundary 
data is the value of the function over part of the boundary. If 
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the boundary data is sufficiently smooth, we can transform the 
homogeneous equation with non -homogeneous boundary conditions to 
a non -homogeneous equation satisfying homogeneous boundary condi- 
tions. The problem can be stated in the following form then. 


F(y)cp xx - ^ = f x (x,y) on fl 

(4.2) 

a ^2. + b ^2. = 0 on dfi 
os on 



where a or b or both may be zero on part of dft. It is 
possible, also, that the stronger condition dcp/Ss = Scp/Sn = 0 
may be imposed on some portion of dfl. 

4.2 Tricomi Equation in Symmetric Positive Form 

It is desired to express (4.2) as a system of first order 
differential equations. Following Friedrichs, we can do this 


by letting u-|_ = cp x , Ug = q> y > and u 

bility condition, cp = q) , we have 

«y ^ 


^F(y) 0\ ( 0 -1 


0 


v-l 


0 



Using the compati- 



(4.3) 


This equation is symmetric, but not positive, since G = 0. To 
make (4.3) symmetric positive, we can multiply by a 2 X 2 matrix, 
B. In order to keep the coefficient matrices of ^u/Sx and 
du/dy symmetric, B must be of the form 

(b cF(y) 
b 


B = 


(4.4) 


where b and c are functions of x and y. 
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Equation (4.3) can now be expressed in symmetric positive 
form by 

Ku = f on ft 

where 

Ku = <x x Ux + (a^uJx + ayhy + (a7u)y + Gu 
^bF(y) cF(y) 

^cF(y) b 

/cF(y) b 

b c 


(4.5) 


x 1 

0U X = 7T 


a y = _i 
a 2 



> (4.6) 


by - c x F(y) 


c - b 

y ^ 


f = 


'bf., 




For the proper choice of functions b and c, G will be positive 
definite, resulting in symmetric positive K. The specific choice 
of b and c depends on the shape of ft and on the boundary 
conditions which are specified. 

It is possible that B may be singular in ft, however, this 
does no harm if B is singular only along a line. 

4.5 Admissible Boundary Conditions 

Let n = (rix,ny) be the outer normal along dft. Then 
6 = n 'a or 



(4.7) 



F(y)(n x b - rye) rycF(y) - ryb N 


n x cF(y) - riyb n x b - rye 


Friedrichs [ i ] noted that the quadratic form u • |3u may be 
written 

-2 „ 2 -p ,, v . \2 _ ( ^,2 _ ».2 ^ ( v ,, a . \2 


(b 2 - c d F(y))(a v .u 1 -ryu 2 ) 2 - (n 2 -F(y)nf)(bu 1 +cu 2 )‘ 

u • 6u — iL ttht — t 

M 2(bn x + cry) 


(4.8) 


From this we can easily specify the boundary matrix, p, so that 
admissible boundary conditions result. Let p be defined so that 

|b 2 -c 2 F(y) Kiyy -ryug) 2 + |ry -F(y)n 2 | (bu-]_ +cu 2 ) 2 ^ ^ 

u M- u ~ 0 i, , i (4= .9) 

2|bn x + cry | 

Of course p is non-negative definite. Also, T)(p - p)©t)(p + p)=R 2 , 
so that the boundary condition Mu = (p - J3) u=0 is admissible. 
Thus we can always obtain admissible boundary conditions. How- 
ever, since b and c can be chosen subject only to the con- 
straint that G is positive definite, we can have a wide variety 
of possible boundary conditions. 

The actual boundary conditions for cp are determined by the 
signs of bry + cry, b 2 - c 2 F(y), and n 2 - F(y)n 2 . For the 
elliptic part of ft, y < 0, we have F(y) < 0; hence both 
b 2 - c 2 F(y) > 0, and n|. - F(y)n x > 0, so that the boundary 
condition is determined solely by the sign of bn + cn . 

t/ 

Suppose that bry + cry '< 0. Then 

|b 2 - c 2 F(y) | (n^ - n^) 2 
u • Mu = u • (p - 0)u = P T *. T" 1 


(4.10) 
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so that Mu = 0 implies n y up = n^ug* In terms of cp this means 
that <P x /9y = r^/ny. ~ -(dy/ds)/(dx/ds), where y = y(s), x = x(s) 
are boundary coordinates as a function of arc length, s, along 
Hence 


<p & + <p $L m o 

ds T x ds ^y ds 


(4.11) 


so' that <p is constant along dft when bn^. + criy < 0 in the 


elliptic region. On the other hand, if br^ + criy > 0, we have 

~2 T1 / \ 2 | /-u,, .1 \Z 

(4.12) 


n • Mu - Ifj ' F (y)P-xl (bu l * cu 2>' 
|bn x + cnyl 


so that Mu = 0 implies that bu^ + cug = 0, so that dcp/dp = 0, 
where p is in some non -tangential direction. Thus, for the 
elliptic region we generally have a single boundary condition 
corresponding to the usual elliptic type of boundary condition. 

In the hyperbolic region the boundary conditions depend on 
whether the magnitude of the boundary slope is greater than, less 
than, or equal to the magnitude of the slope of the characteristic 
curve. For equation (4.2), the characteristics satisfy the differ- 
ential equation 


%L 

dx 



VHy) 


(4.13) 


Suppose that the boundary is tangent to a characteristic, then 


fx _ _ d£ 

■V ’ 
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so that 

= F(y)n| (4.14) 

Suppose that a portion of 60 is a left running characteristic, 
so that ry = V F (y) “x* ^en, from (4.8) 

u * pu - — * c V F nT)(b - c-v/F(y))(n y u 1 - n x u 2 ) 2 

2(b + c-y/F(y) )n x 

2 

„ (h - cVF(y))(n y u 1 - n x u 2 ) 

2n x 

Suppose that this portion of 60 is a right boundary, so that 
n^. > 0, then, if b < c^F(y), equation (4.1l) holds and 
dcp/ds = 0, but if b > c-/F(y ) , p = 8 and no boundary condition 
is imposed. Similarly, along a right running characteristic, 
the same types of boundary conditions are determined by the sign 
of b + cVF(y) . 

2 2 

If the boundary is not characteristic, we have riy > F(y)n£ 
if the magnitude of the boundary slope is less than the magnitude 
of the characteristic slope, and vice versa. Thus, the particular 
boundary condition is determined by the signs of all three terms, 
n 2 - F(y)n^, b 2 - c 2 F(y), and br^ + cn y . 

The various boundary conditions implied by a choice of b 
and c are summarized in Table I. It can be seen that there 
is considerable choice in the type of boundary which can be 
specified by the proper choice of the functions b and c. 
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TABLE I. - SUMMARY OF BOUNDARY CONDITIONS 


Elliptic part of D (y < 0) 

Boundary condition 

Condition on b and c 

P = o 

ds 

bn x + cny < 0 

o 

ti 

#ia 

bn x + cny > 0 


Hyperbolic part of £2(y > 0) 


Boundary condition 

Type of boundary 

Conditions 

on 

b and c 

d9 _ 
ds 

0 

n y = ^F(y)n x ,n x >0 

b ^ cVF(y) 



none 


n y * V F (y) n x^ n x > o 

b > c-/F(y) 



none 


ny= -n/FfyV^ny.X 0 

bb< - cVRF) 


d<p _ 
ds 

0 

n y = -V F (y) n x^ n x < 0 

b > - c-y/F(y) 


dcp _ 
dp " 

0 

ny > F(y)n x 

b 2 >c 2 F(y) 

and 

bn x +cny>0 

dqp _ 
ds 

0 

n| > F(y)'n 2 

b 2 > c 2 F(y) 

and 

bn x + cn y < 0 

d(p _ 
ds ~ 

dn 

n| > F(y)n 2 

b 2 < c 2 F(y) 

and 

bn x + cny> 0 

none 


n| > F(y)n x 

b 2 <c 2 F(y) 

and 

bn x + cny < 0 

none 


n y < F(y)n x 

b 2 > c 2 F(y) 

and 

bn x + cny > 0 

dcp = 

12 = 0 

n| < F(y)n 2 

b 2 > c 2 F(y) 

and 

bn x + cny < 0 

ds 

dn 




dcp _ 
ds 

0 

n y < F(y)n x 

b 2 < c 2 F(y) 

and 

bn x + cn y > 0 

dcp 

0 

n 2 < F(y)n x 

b 2 < c 2 F(y) 

and 

bn x + cny < 0 


Note: p denotes the distance in some non-tangential direction. 
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4.4 Sample Problem 

A simple choice of b 
positive definite in ft, if 

b = -b 

c = c Q , where 


and c which will result in 
F* (y) > 0, is 
0 - b-,x, b x > 0 

bnF(y) 

c n > - - L. ILL in ft 
0 F' y 


G 


being 


(4.16) 


Then 


A>iF(y) + c 0 F' (y) 

\ ° 



(4.17) 


which is obviously positive definite. 

To show the type of boundary conditions which may result, 
consider the case F(y) = y, so that 


G 



'V + c 0 

0 



(4.18) 


The characteristics in this case satisfy one of the equations 


dx 


± 


1 

VSF 


which can be solved to obtain the characteristic equation, 


y 3 = I (x - x 0 ) 2 


(4.19) 


where xq is the point on the x-axis intersected by the 
characteristic . 

As ' an illustration, suppose that ft is the region shown in 
figure 2, which is bounded by two characteristics in the hyper- 
bolic region and by a curve satisfying riy <[(b Q + b 1 x)/c Q ] n x 
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in the elliptic region. It is assumed that bo/bp is chosen 

p p 

large enough so that the parabola (bpx + bp) = c^y lies entirely 
to the left of fl, as indicated in figure 2. The boundary condi- 
tion is dcp/ds = 0 for the elliptic portion, 

Tp, of dft, and for one characteristic, Pg, with no boundary condi- 
tion on the other characteristic, T3. This is known as a Tricomi 
problem. Variations are possible with Tg and replaced by 

several characteristics. This type of problem is discussed by 
Bers [ 11 ] , p . 88 . 


It is worthwhile noting that the solution obtained by the 
finite difference solution of the symmetric positive form of the 
Tricomi equation consists of derivatives of the stream function, 
which corresponds to velocities in the physical problem. Hence, 
even though we have a convergence rate which is less than CKh 1 / 2 ), 
it is essentially equivalent to a convergence rate of 0(h 3 / 2 ) 
if the original second order equation were solved directly for 


the stream function. 


CHAPTER V 


A NUMERICAL EXAMPLE 
5.1 Description of Problem 

A numerical solution to a Tricomi equation was calculated 
using the finite difference scheme of Chapter II. The accuracy 
of the solution was checked by using a problem for which an 
analytical solution is known. 

The Tricomi equation can be put in symmetric positive form 
as indicated in the last chapter, as given by equations (4.5) 
and (4.6). The region A chosen is indicated in figure 3. 

The choice of b and c are 

b = -3 - x 

\ (5.1) 

c = 2 J 

which gives b x = -1, and b y = c x = c y = 0. We choose F(y) = y. 
Using this in (4.6) we have 



which is positive definite in A. 

We now check to see what the admissible boundary conditions 
are from Table I. For the hyperbolic part of A(y > 0) , we need 
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2 2 * 2 

to know the sign of b p c y = (x + 3') - 4y. From figure 3 we 
see that b^ > c^y in ft, hence, since b < 0 in ft, we have 
± cVy < -b in ft. From Table I, we have dcp/ds = 0 on Tp 
and no boundary conditions on Tg. For the elliptic part of 
ft(y < 0) we need to check the sign of bn^. + cny. Along T g 
we have = 1, riy = 0, so that bn^. + cny = -(x + 3) < 0 
along Tg • Hence, the admissible boundary condition along Fg 
is dcp/ds = 0. Next we check r 3 . Then n x = 0 and ny = -1, 
and bn x + cn y = -2 < 0, so that dq>/ds = 0 along r 3 . Finally, 
along r 4 , n^. = -1, riy = 0, giving bn x + cny = x + 3 = 2 > 0, 
since x = -1. Hence dcp/dp = 0 along r 4 , where p is some 
non-tangential direction. To find the specific direction, we go 
back to equation (4.12) which holds in this case. We see that 
Mu = 0 implies that bup + cug = -2cp x + 2 cp y =0 or qp x = <p y . 

Hence p is in a direction sloping downward at 45°. We summarize 
the boundary conditions : 



Boundary Condition 

r l 

sit 

II 

o 

r 2 

o 

II 

SI-3 

r 3 

o 

II 

r 4 


r 5 

None 


A simple, but non-trivial function satisfying these boundary 
conditions is easily obtained by choosing a function which is 
zero along Tp, Tg, fg, and r 4 , with the normal derivative also 
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zero along Pq. These requirements are met by 

cp(x,y) = (x + l) 2 (x - l) (y + l)(4y 3 - 9(x - l) 2 ) (5.3) 

The function fq is determined by calculating ycp^ - <p = 
which gives 

fl(x,y)» y(y + l)[(4y 3 -9(x -l) 2 (6x+2) -18(x 2 -l)(7x -l)] 

- 24(x + l) 2 (x- l)y(2y + l) (5.4) 


The functions for which we are solving are then 
9 X = (x +1) (y + l)[ (4y 3 - 9(x - l) 2 ) (3x - 1) -18(x + l) (x - 


l) 2 ] 


cpy — (x +l) 2 (x - l) [l6y 3 + 12y 2 - 9(x -l) 2 ] 


and from (4.6) we have 



(5.5) 


(5.6) 


(5.7) 


We need to evaluate the matrix p along all boundaries , with 
p defined by equation (4.9) . A straightforward calculation gives 
the following values for p . 
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This gives the information necessary to calculate the coef- 
ficients of the finite difference equation, which is 

2 L j,k P j,k u k + 2 + A j G j u j = A j f j ^ 5,< 

Equation (5.8) holds for every mesh point, x i? in the set of mesh 
points. For simplicity a uniform mesh was used, as indicated in 
figure 4. It will be noted that mesh points outside of ft were 
used. A solution was calculated for two different mesh spacings, 
h = 0.2 and h = 0.1. The finite difference equation was solved 
in each case by the block tridiagonal algorithm mentioned in 
section 2.5. Since the analytical solution, u, is given by (5.5) 
we can calculate ||u h - as well as the maximum value over 

all mesh points of the maximum component of the error. 





Figure 4. - Mesh point arrangement for numerical example. 
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5.2 Description of Numerical Results 

Theorem 2.1 assures us of essentially Oth 1 / 2 ) convergence in 
2 

the L norm. Unfortunately this does not assure us of point- 
wise convergence. As indicated in the proof of Theorem 2.1, the 
finite difference equations can be expected to be less accurate 
when the polygons, Pj, are not uniform rectangles. This was the 
case in the numerical example. The result was poor accuracy near 
the hyperbolic boundary segments, and Ig. In going from 
the coarse mesh (h = 0.2) to the fine mesh (h = 0.1), the L 2 
error was reduced from 6.06 to 5.30 which is not unreasonable 
with the Odi 1 / 2 ) convergence rate. However, the maximum error 
actually increased from 33.5 to 60.9 indicating pointwise diver- 
gence. The horizontal line (y = 0.75) along which the finite 
difference solution for the finer mesh has the poorest agreement 
with the analytical solution is plotted in figure 5. It is seen 
that the finite difference solution has large oscillations with 
a "wild" point at the end of the line. 

2 

All this is not quite as bad as it seems, though, since L 
convergence with pointwise divergence means that the divergent 
points will occur as sharp peaks. Therefore, it' can be expected 
that a smoothing operation would give great improvement inthe 
results. With this in mind, a simple smoothing procedure was 
tried. Since most smoothing methods are for one -dimensional 
functions, the solution was smoothed by lines, first along vertical 
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lines and then along horizontal lines. The method of smoothing 
used is similar to a method suggested by Hamming, p. 314, [12]. 
If it is desired to smooth equally spaced data, , we can 
define the smoothed data, k=i> by 


y o " 


y 0 + 2y l - y 2 


= y k-l + 2y k + y k+l .... „ „ 

y k — 4 : ~ * for k — 2,3, . . . , n - 1 


- . - y n-2 + 2y n-l + y n 
2 


The result of applying this smoothing procedure to the 
solution based on the finer grid (h = O.l) was to reduce the L 2 
error from 5.30 to 2.07. The maximum error was reduced from 
60.9 to 13.8. This maximum error was at a point lying outside of fl, 
the maximum error for a mesh point in ft was 6.4. The improve- 
ment obtained by this smoothing procedure is indicated by figure 6, 
which shows the horizontal line with poorest agreement after 
smoothing. The solution after smoothing along a more typical 
horizontal line is shown in figure 7. 

It s-hould be emphasized that the smoothing procedure used here 
was very simple and that most likely better results could be 
obtained with other smoothing methods. For example, Lanczos [13], 
gives several smoothing methods, both local and global (through the 
use of truncated Fourier series) . 
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